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Preface 


This  report  describes  the  work  performed  under  the  U.S.  Army  Engineer 
Waterways  Experiment  Station  (WES)  Laboratoiy  Discretionary  Research  Program 
titled  "Fidl  Waveform  Inverse  Modeling  of  Ground  Penetrating  Radar  Data."  The 
research  was  performed  by  Drs.  Janet  E.  Simms  and  Dwain  K.  Butler,  Earthquake 
Engineering  and  Geosciences  Division  (EEGD),  Geotechnical  Laboratory  (GL). 
Outside  technical  assistance  was  provided  by  Mr.  Michael  H.  Powers,  United  States 
Geological  Service,  Denver,  CO.  Ms.  Lori  Davis  assited  in  preparing  this  report. 
The  work  was  performed  imder  the  direct  supervision  of  Mr.  Joseph  R  Curro,  Jr., 
ChieC  Engineering  Geophysics  Branch,  EEGD.  The  work  was  performed  under  the 
general  supervision  of  ^s.  A.  G.  Franklin,  Chief,  EEGD,  and  William  F.  Marcuson 
in.  Director,  GL. 

At  the  time  of  publication  of  this  report.  Director  of  WES  was  Dr.  Robert  W. 
Whalin.  Commander  was  COL  Bruce  K.  Howard,  EN. 


The  contents  of this  report  are  not  to  be  used  for  advertising,  publication, 
or  promotional  purposes.  Citation  of  trade  names  does  not  constitute  an 
official  endorsement  or  approval  of  the  use  of  such  commercial  products. 


Conversion  Factors, 
Non-SI  to  Si  Units  of 
Measurement 


Non-SI  units  of measurement  used  in  this  report  can  he  converted  to  SI  units  as  follows: 


1  Introduction 


Background  and  Applications 

Ground  penetrating  radar  (GPR)  is  an  electromagnetic  (EM)  geophysical  method 
that  has  emCTged  as  a  versatile  high  resolution  tool  for  a  variety  of  applications. 
When  GPR  systems  became  commercially  available  in  the  1970's,  the  capabilities  of 
GPR  were  easUy  oversold  to  a  market  eager  for  “high-tech”  solutions  to  difficult 
subsurface  problems.  In  certain  geologic  settings,  even  the  veiy  early  GPR  systems 
could  produce  real  time  graphic  records  (i.e.,  no  processing  required)  which 
resembled  snapshots  of  the  subsurface.  This  tendency  to  oversell  the  method  led  to 
frequent  application  of  GPR  to  inappropriate  sites  and  objectives.  Another  problem 
was  the  use  of  GPR  by  inexperien^  personnel.  Fortunately,  the  capabilities  and 
limitations  of  GPR  are  now  better  understood  by  geophysicists  and  the  geotechnical 
and  geoscience  community  as  a  whole,  and  the  method  is  being  more  rationally 
exploited.  In  cases  where  GPR  has  a  sufficient  depth  of  investigation  capability  to 
address  a  specific  problem,  it  has  the  greatest  horizontal  and  vertical  resolution  of 
any  geophysical  method  that  can  be  applied  to  the  problem. 

GPR  is  applicable  to  a  broad  spectrum  of  interests  and  problems  in  geotechnical 
engineering  and  geoscience  (Butler  1992).  Some  of  the  areas  of  application  include 
the  following:  geological  mapping  (Butler  et  al.  1990;  Davis  and  Aiman  1989), 
foundation  investigations  (Butler  et  al.  1989),  cavity  and  sinkhole  mapping  in  karst 
areas  (Ballard  1983),  ground  water  investigations  (Bohling  et  al.  1989), 
archaeological  studies  (Butler  et  al.  1994),  detection  and  mapping  of  buried  cultural 
features,  e.g.,  waste  disposal  trenches,  utilities,  barrels,  underground  storage  tanks, 
unexploded  ordnance,  etc.  (Sharp  et  al.  1990;  Llopis  and  Simms  1993),  contaminant 
detection  and  mapping  in  the  subsurface  (Brewster  and  Annan  1994),  and  structural 
studies  of  buildings,  bridges,  dams,  etc.  In  a  military  context,  key  areas  of 
application  are  pavement  assessment  and  evaluation  (Maser  and  Scullion  1992), 
airfield  evaluation,  tunnel  detection,  concrete  structure  investigations,  installation 
restoration,  and  foundation  investigations. 

Most  commercially  available  GPR’s  are  short  pulse,  time  domain,  EM  wave 
propagation  systems.  The  GPR  systems  use  closely  spaced  EM  transmitters  (Tx) 
and  receivers  (Rx),  and  the  common  graphic  data  display  of  GPR  waveforms  side- 
by-side  resembles  a  common  depth  point  (CDP)  seismic  section  (Figure  1).  Jn  the 
cartoon  in  Figure  la,  the  Tx  and  Rx  are  both  contained  in  the  box  indicated  as 
“antennae”;  for  many  GPR  systems  the  Tx  and  Rx  are  separate,  and  the  spacing 
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between  them  can  be  varied.  When  interpreting  the  GPR  record,  however,  factors 
must  be  considered  which  are  not  relevant  in  seismic  reflection:  scattering  from 
gravel,  cobbles,  microcrack  clusters,  and  other  small  geologic  inhomogeneities; 
backscatter  signatures  from  overhead  powerlines  and  surface  structures  which  can 
masquerade  as  subsurface  features  in  the  GPR  record;  and  minor  changes  in  the 
dielectric  properties  of  soil,  rock  or  water  which  can  appear  as  major  GPR 
signatures  but  do  not  affect  other  physical  properties,  such  as  seismic  velocity. 
Advantages  of  GPR  are  high  vertical  and  horizontal  resolution,  rapid  survey 
capability,  and  near-real  time  data  interpretation  in  many  cases.  Major 
disadvantages  are  the  site  specific  applicability  and  limited  depth  of  investigation 
capability  (generally  less  than  15  m). 


Justification  and  Reievance 


Ground  penetrating  radar  surveys  are  often  interpreted  in  a  qualitative  manner  to 
identify  anomalies.  This  method  of  interpretation  is  used  to  determine  the  existence 
of  a  feature,  when  specific  properties  of  tee  feature  are  of  no  interest.  For  example, 
in  tee  case  of  tee  hypothetical  GPR  cross-section  in  Figure  la,  tee  extent  of  tee 
interpretation  typically  is  limited  to  (1)  identification  of  radar  events  corresponding 
to  subsurface  reflectors  (interfaces)  and  (2)  determination  of  depths  to  tee  interfaces 
assuming  an  average  EM  wave  propagation  velocity.  However,  GPR  data  can  also 
be  interpreted  quantitatively  to  give  more  definitive  depths,  geometry,  and  EM 
(physical)  properties  of  tee  subsurface  material.  Existing  analysis  procedures  are 
simplistic,  making  use  of  only  a  small  amount  of  tee  information  contained  in  tee 
received  GPR  waveforms,  and  do  not  consider  electromagnetic  wave  attenuation, 
dispersion,  scattering,  geometric  spreading,  and  interface  refiection/transmission 
losses.  The  aim  of  this  research  is  to  develop  interpretation  procedures  which  utilize 
tee  full  GPR  waveform  (Figure  lb)  for  determining  layer  thicknesses  and  tee  EM 
properties  of  subsurface  materials. 


Scope  of  Work 

This  woric  consisted  of  three  phases:  1)  review  existing  capability  for  full 
waveform  forward  modeling  of  GPR  signals,  2)  develop  a  catalog  of  theoretical 
GPR  receiver  waveforms  for  representative  layered  models  (pavements,  concrete 
with  air-  and  water-filled  cavities,  etc.)  of  interest  using  typical  physical  property 
values,  and  3)  examine  tee  feasibility  of  developing  an  inverse  modeling  procedure 
to  determine  physical  properties  by  comparing  full  waveform  forward  solutions  to 
measured  GPR  signatures. 

The  first  phase  involved  identifying  available  GPR  modeling  programs  and 
determining  which,  if  any,  are  suitable  for  performing  full  waveform  forward 
modeling  of  GPR  signals  for  use  in  an  inverse  modeling  program  or  procedure.  The 
primary  criteria  are  teat  tee  model  assumptions  must  be  realistic  and  defensible  and 
teat  tee  major  EM  phenomena  (spreading,  attenuation,  dispersion, 
reflection/transmission,  scattering,  multiples)  be  included.  One  GPR  forward 
modeling  program  met  tee  requirements,  ''GPRMODv2''  (Powers  and  Olhoefl 
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1995).  Since  this  project  began,  GPRMODv2  underwent  a  major  revision  and  the 
authors  were  cooperative  in  incorporating  special  program  features  to  facilitate  this 
research.  The  concepts  of  GPR  surveying  and  the  key  theoretical  aspects  of  EM 
wave  propagation  in  a  layered  media  are  developed  in  Chapter  2.  The  major 
assumptions  and  conditions  of  the  forward  modeling  program  will  be  examined  in 
Chapter  3. 

For  Phase  2  of  this  project,  the  forward  modeling  program  GPRMODv2  was 
tised  to  generate  models  representative  of  pavement  structure  and  concrete  with  air- 
or  water-filled  cavities.  Parameters  were  varied  for  both  two  and  three  layer  models 
to  investigate  their  affect  on  the  resulting  GPR  waveform. 

Finally,  for  Phase  3,  the  requirements  of  developing  an  inverse  GPR  modeling 
capability  are  surveyed.  First  the  procedures  for  manually  adjusting  the  parameters 
of  a  subsurface  model  in  an  iterative  manner  to  match  model  output  with  field  GPR 
records  is  considered.  Next  the  possibilities  and  advantages  of  simplifying  the 
invasion  process  by  constraining  and/or  reducing  the  number  of  model  parameters 
is  addressed.  The  number  of  model  parameters  can  be  reduced  by  simply 
eliminating  some  parameters  fi'om  consideration  in  the  inversion  or  by  determining 
some  parameters  in  advance  by  other  means.  If  some  parameters  are  known  or 
assumed  for  the  subsurface  model,  then  they  can  be  fixed  or  constrained  during  the 
invCTsion,  thereby  simplifying  the  inversion  process.  Parameters  can  be 
independently  determined  by  laboratory  testing  on  soil/rock  samples  in  the 
laboratory,  from  drilling  logs,  fi'om  the  results  of  electrical  resistivity  and 
electromagnetic  soundings,  fi'om  wide-angle  GPR  reflection/refiaction  surveys,  and 
fi'om  push  probe  time  domain  reflectometiy  (TDR).  The  possibility  of  cross- 
correlation  of  the  GPR  waveform  with  the  transmitter  waveform  or  deconvolution  of 
the  GPR  waveform  to  determine  number  and  location  of  interfaces  as  the  first  phase 
of  inversion  is  studied.  Finally,  the  type  of  inversion  procedure  most  likely  to 
succeed  will  be  determined;  global  inverse  solution,  such  as  typically  used  for 
electrical  resistivity  and  fi'equency-  and  time-domain  EM  inversion,  or  a  "layer¬ 
stripping"  inversion  process. 
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2  Electromagnetic  Wave 
Theory 


Ground  penetrating  radar  is  an  electromagnetic  (EM)  method  which  operates  in 
the  megahertz  (MHz)  or  gigahertz  (GHz)  frequency  range.  At  these  high 
frequencies  displacement  currents  (currents  associated  with  charges  which  are 
constrained  from  moving  continuously)  dominate  and  energy  propagates  into  the 
ground  as  a  wave.  An  electromagnetic  pulse  transmitted  into  the  earth  imdergoes 
refraction,  reflection,  scattering,  and  dispersion.  Contrast  in  the  dielectric 
permittivity  at  layer  boundaries  causes  Ihe  EM  wave  to  be  reflected  and  refracted. 
The  dielectric  permittivity  is  the  proportionality  factor  relating  the  displacement 
current  to  energy.  Since  electromagnetic  fields  consist  of  both  electric  and  magnetic 
fields,  any  property  of  the  geologic  material  which  affects  either  of  these  fields  will 
also  affect  the  propagation  of  the  EM  wave  in  the  subsurface.  Generally,  the  elec¬ 
trical  properties  of  the  soil  and  rock  have  a  greater  variation  than  the  magnetic 
properties  and,  therefore,  a  greater  influence  on  the  EM  wave  propagation.  Soil 
conductivity,  which  is  usually  controlled  by  the  moisture  content,  is  a  major  factor  in 
determining  if  GPR  can  be  used  successfully  at  a  site.  High  conductivity  soils,  such 
as  those  with  high  clay  content,  can  significantly  attenuate  the  EM  signal  and  render 
GPR  virtually  useless. 

Propagation  Velocity,  Attenuation,  and  Wavelength 

The  governing  wave  equation  appropriate  for  ground  penetrating  radar  is  derived 
from  Maxwell’s  equations  (Balanis  1989,  chapters  1  and  3).  For  a  plane  wave 
moving  through  a  source-free  material  it  is 

V2E  =  -k2E  or  V2E  =  y^E, 

where  k  is  the  complex  wave  vector,  y  is  known  as  the  propagation  vector,  E  is  the 
electric  field  vector,  and  is  the  Laplacian  operator.  Both  k  and  y  are  used  in  the 
literature,  and  the  relationship  between  them  is  simply  y=ik  where  1=^ - 1 .  The 
material  can  be  dispersive  (velocity  varying  with  frequency)  and  lossy, 
inhomogeneous,  anisotropic,  or  none  of  these,  depending  on  how  the  wave  vector 
(or  number)  or  propagation  vector  (or  constant)  is  characterized.  In  the  most 
general  terms,  a  solution  to  this  equation,  which  describes  the  plane  wave  moving 
through  the  material,  is 
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E(r)-e.E,e*'»''> 

where  r  is  the  position  vector,  e,  is  a  unit  vector  in  the  direction  of  the  electric  field, 
and  Eq  is  an  initial  amplitude  of  the  electric  field.  If  a  harmonic  time-dependence  is 
assumed,  as  in  a  propagating  EM  wave,  this  general  solution  becomes 

E(r)  =  e.Real  [e‘“‘Eoe'*' 

If  the  material  is  assumed  to  be  homogeneous,  then  its  properties  do  not  vary  as 
a  function  of  position.  Almost  all  rocks  and  soils  have  some  heterogeneities,  but  on 
a  scale  similar  to  the  GPR  wavelength,  they  can  generally  be  divided  into  volumes 
of  homogeneous  materials  with  the  property  differences  occurring  at  defined 
boundaries.  If  the  material  is  isotropic,  then  the  properties  are  not  direction 
dependent.  Anisotropy  is  not  imcommon  in  rocks,  as  most  crystalline  structures 
show  distinct  anisotropy,  but  again  on  the  scale  of  most  GPR  wavelengths  an 
assumption  of  isotropy  is  not  unrealistic.  The  material  can  also  be  assumed  linear, 
meaning  its  properties  do  not  vary  with  the  amplitude  of  the  electric  or  magnetic 
field.  This  is  true  at  the  low  magnitudes  of  the  fields  produced  by  ground 
penetrating  radar  systems.  See  Keller  (1987)  or  Balanis  (1989,  chapter  2)  for  more 
information  on  these  material  property  assumptions. 

As  an  EM  wave  moves  through  a  material,  a  given  location  experiences  an 
electric  and  magnetic  field  that  reverses  polarity  with  time.  When  the  material  is  in 
the  presence  of  an  electric  field,  some  energy  is  lost  as  charge  transport  takes  place 
according  to  the  conductivity.  Some  energy  is  stored  by  charge  polarization 
processes  according  to  the  measure  of  dielectric  permittivity.  However,  different 
polarization  processes  take  different  amounts  of  time  to  complete,  and  the  frequency 
of  the  EM  wave  will  determine  how  much  energy  is  stored  for  a  given  polarity  cycle. 
This  results  in  a  fi'equency  dependent  dielectric  permittivity.  Each  polarization 
process  also  involves  some  amount  of  energy  loss  through  brief  conductive  losses 
while  the  charges  are  moving  to  their  polarized  state.  This  is  the  displacement 
current  loss,  sometimes  called  an  alternating  conductivity,  or  a  dielectric  relaxation 
loss,  and  can  be  characterized  by  making  the  dielectric  permittivity  complex 
(Olhoeft  1981).  The  magnetic  permeability  of  most  rocks  and  soils  is  not  much 
different  fi'om  the  magnetic  permeability  of  fi^ee  space,  but  when  iron  is  present  it 
can  also  be  fi'equency  dependent  and  complex  with  magnetic  relaxation  effects 
(Olhoeft  and  Strangway  1974;  Olhoeft  and  Capron  1994). 

For  isotropic,  homogeneous  materials,  the  wave  vector  or  propagation  vector 
reduce  to  complex,  fi'equency  dependent  values  called  the  wavenumber  or 
propagation  constant.  Depending  on  the  properties  of  the  material,  the  wavenumber 
and  propagation  constant  may  be  further  reduced.  In  free  space,  where  the 
conductivity  is  zero  and  propagation  exists  with  no  material  losses,  the  wavenumber 
is 


where  o)  is  the  angular  fi'equency,  po  is  the  magnetic  permeability  of  fi-ee  space 
(47tx  10'^  Henries/meter),  and  Cq  is  the  dielectric  permittivity  of  fi'ee  space 
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(8.854x10’*^  Farads/meter).  The  solution  for  a  plane  wave  traveling  in  the  z 
direction  reduces  to 

E(z)  =  e^Eo  cos[a)(t-z/c)] 

where  the  propagation  velocity  is 

V  =  ^  -  =  c  =  3  X 1 0*  meters/second, 

^0 

which  is  the  speed  of  light  in  free  space. 

It  is  often  necessary  for  simplification  to  assume  that  the  rocks  and  soils  being 
investigated  with  GPR  are  homogeneous,  isotropic,  linear,  and  source-free. 
However,  they  may  be  dispersive  and  lossy,  which  can  be  handled  without  too  much 
complexity.  The  wavenumber  can  be  written  as 


(  Mr"  /’Mr  )  K  )  -  Mo  (  Mr  /'Mr  )  % 


where  Pj.  and  Ej.  are  the  complex,  frequency  dependent  relative  permeability  and 
permittivity,  and  Og  is  the  static  (zero  frequency)  conductivity.  Relative  permittivity 
and  permeability  are  ratios  of  the  real  property  value  with  respect  to  the  value  of 
free  space  as  in 


t  and  =  = 

Mo 

where  e/,  e‘,  p,’,  and  p^’  are  the  real  and  imaginary  parts  of  the  relative,  complex 
permittivity  and  permeability,  respectively.  For  this  case,  the  wavenumber,  k,  is  a 
complex,  frequency  dependent  value.  It  can  be  expressed  in  terms  of  its  real  and 
imaginary  parts  as 


)fc(ci))=p(G))-ra(a)) 

where  P  is  the  real  part  of  the  wavenumber  and  a  is  the  imaginary  part.  Both  of 
these  terms  vary  with  frequency.  In  these  terms,  the  propagation  constant  is  simply 

Y(o>)=oc(ci>)+;P(a)) . 

At  this  point  the  terminology  choice  between  ^  and  y  is  irrelevant,  as  P  is  the  phase 
parameter  and  a  is  the  attenuation  parameter  for  either  choice.  The  wave  equation 
solution  for  a  plane  wave  moving  through  this  material  in  the  z  direction  is 

E(z)  =  e,(Eo  e"®^  cos(ot  ±  Pz)  .  .  ' 

This  wave  is  characterized  by  a  velocity  of 
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a  wavelength  of 

1  271 

and  an  amplitude  change  in  decibels  per  meter  of 
a((o)=  -8.686  a((i)) 

due  to  conductive  and  relaxation  losses.  The  frequency  dependent  parameters  a  and 
P  can  be  expressed  in  terms  of  the  material  properties  as 


\^ere 


A 


and 


A  convenient  way  to  handle  both  the  frequency  dependent  and  complex  nature  of 
the  permittivity  and  permeability  is  through  the  Cole-Cole  expression  (Cole  and 
Cole  1941;  Keller  1987;  Olhoeft  and  Capron  1994).  The  formulas  are 


(e/-ej 

(1  +(70)t/‘) 


and 


11-  + 


(1 


where  e_  and  e,  and  p,  are  the  high  and  low  frequency  limits  of  permittivity  and 
permeability,  respectively;  x*  and  are  time  constants  of  relaxation,  and  and 
are  Cole-Cole  time  constant  distribution  parameters  (Cole  and  Cole  1941).  In  each 
layer,  for  both  permittivity  and  permeability,  four  values  are  given  (low  frequency 
limit  high  frequency  limit,  time  constant,  and  distribution  parameter).  If  the  high 
and  low  limits  are  equal  to  each  other,  the  property  value  is  real  and  independent  of 
frequency  (and  the  time  constant  and  distribution  parameter  are  irrelevant).  When 
the  high  and  low  frequency  limits  are  different,  the  functions  behave  as  shown  in 
Figure  2.  The  time  constant  of  relaxation,  Xg  or  x^,,  controls  where  in  frequency  the 
decrease  in  the  real  part  steps  from  e,  to  e„,  which  is  also  the  frequency  where  the 
imaginary  part  is  at  a  maximum.  The  slopes  on  the  step  down  in  the  real  part  and 
the  flanks  of  the  imaginary  part  are  controlled  by  the  distribution  parameter,  or 
a^,  which  varies  from  1.0  to  0.0  for  steepest  to  smoothest  slope.  When  either  or 
Cjj  is  zero  the  corresponding  permittivity  or  permeability  is  a  constant  value  midway 
between  the  high  and  low  frequency  limits. 

Using  the  Cole-Cole  functions  means  ten  parameters  are  required  to  describe 
each  earth  material,  four  each  for  the  complex  magnetic  permeability  and  the 
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ccMnplex  dielectric  permittivity,  plus  a  static  conductivity  and  a  thickness.  In  the 
frequmcty  range  from  above  100  MHz  to  less  than  1  GHz,  soils  and  rocks  often 
have  very  small  relaxation  losses  (Olhoeft  and  Capron  1994).  By  disregarding  the 
complex  nature  of  the  material  properties  and  assuming  that  no  iron-bearing 
matmals  are  present,  which  are  reasonable  assumptions  for  many  cases,  the  number 
of  parameters  can  be  reduced  to  three  per  material.  They  are  a  real,  relative 
dielectric  permittivity,  a  real,  static  conductivity,  and  a  thickness.  The  imaginary 
parts  of  the  permittivity  and  permeability  are  assumed  to  be  zero,  and  the  real 
relative  magnetic  permeability  is  assumed  to  be  one.  For  this  case,  the  wavenumber 
is 


and  its  real  and  imaginary  parts  become 


Note  that  because  a  and  P  are  still  ftinctions  of  frequency,  the  wave  traveling 
through  this  material  is  still  dispersive,  but  only  through  conductive  losses. 
However,  when  o/co  Eq  is  on  the  order  of  one  or  less,  which  it  will  be  for 
conductivities  less  than  5  millisiemen/meter  (mS/m)  at  100  MHz  and  less  than  55 
mS/m  at  1  GHz,  then  a  reasonable  approximation  for  p  is 

p-f- 

Because  the  velocity  of  the  wave  is  (o/p  ,  this  results  in  a  frequency  independent 
velocity  of 


which  is  commonly  used  in  ground  penetrating  radar.  The  wavelength  according  to 
this  approximation  is  simply 


where  /is  the  frequency  in  cycles  per  second.  An  approximation  that  can  be  made 
for  a  ^on  Hippel  1954,  page  28)  where 
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results  in  a  formula  for  a  change  in  amplitude  in  decibels  per  meter  according  to: 

a(o))=  -1637— 

fr 


The  vectors  E,  J,  H,  and  B  are  the  electric  field,  the  electric  current  density 
(including  displacement  currents),  the  magnetic  field,  and  the  magnetic  induction, 
respectively.  The  loss  tangent  is  commonly  used  as  a  measure  of  the  loss 
medianisms  in  a  material  with  regard  to  EM  wave  propagation.  However,  separate 
loss  tangents  are  often  defined  for  different  combinations  of  loss  mechanisms.  The 
total  electromagnetic  loss  tangent  is  a/p,  or  the  ratio  of  the  imaginary  to  the  real 
part  of  the  wavenumber.  Physically,  this  is  the  phase  difference  between  the  vectors 
E  and  H  (Lorrain  and  Corson  1970).  The  total  electromagnetic  loss  tangent  is 


( 

tan6^=tan 

em 

\ 


2  J 


a 

p 


The  total  electrical  loss  tangent  includes  conductive  and  dielectric  relaxation  losses, 
but  not  magnetic  relaxation  losses,  and  is  defined  as 


tan  5, 

e 


+ 


0 

(OCo 


e 


r 


Physically,  this  is  the  phase  difference  between  E  and  J.  The  magnetic  loss  tangent 
is 


and  represents  the  phase  angle  between  B  and  H. 


Reflection  and  Transmission  Boundaries 


When  an  EM  wave  encounters  an  interface  between  two  media,  the  contrast  in 
mateial  properties,  the  radius  of  curvature  of  the  interface,  and  the  incident  angle, 
polarization,  and  radius  of  curvature  of  the  arriving  energy  will  determine  the 
amplitudes,  directions,  and  polarizations  of  the  reflected  and  transmitted  energy. 
This  encounter  can  be  very  complex  without  some  simplifying  assumptions. 
Generally,  these  assumptions  include  making  the  incident  wavefi'ont  and  the 
material  boundary  be  planar  surfaces,  and  assuming  only  imiform  plane  waves.  The 
polarization  effects  can  be  eliminated  by  considering  the  problem  in  only  two 
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dimensions,  v^ere  the  incident  wave  is  linearly  polarized  along  the  third  dimension. 
This  is  called  a  horizontal,  perpendicular  or  E  polarization  (Balanis  1989,  chapter 

5). 


Figure  3  shows  a  planar  boimdary  between  two  media,  and  the  rays  associated 
with  the  planar  wavefronts  for  the  incident,  reflected,  and  transmitted  energy.  The 
electric  field  is  oriented  in  and  out  of  the  page  in  this  figure.  For  a  specular  reflector, 
vrii^e  the  surface  roughness  consists  of  height  variations  that  are  small  compared  to 
the  wavelength,  the  angle  of  reflection  is  equal  to  the  angle  of  incidence: 

0,=e, 

The  transmitted  angle  can  be  calculated  from  the  incident  angle  and  the  material 
properties  through  Snell’s  Law  (Balanis  1989;  Lorrain  and  Corson  1970;  Ulaby, 
Moore  and  Fung  198 1).  The  formula  is: 


/ 

sm0,= 


sin0. 


^riiere  and  pi  are  the  wavenumber  and  magnetic  permeability  of  the  incident  layer, 
and  k2  and  p2  the  wavenumber  and  magnetic  permeability  of  the  transmitted 
layw. 

When  the  wavenumber  of  either  material  is  complex,  the  transmitted  angle 
becomes  complex.  Physically,  this  means  for  the  transmitted  wave  the  planes  of 
constant  phase  (direction  of  P,  the  real  part  of  the  wavenumber)  and  the  planes  of 
constant  amplitude  (direction  of  a,  the  imaginary  part  of  the  wavenumber)  are  not 
parallel  (see  Balanis  1989,  page  214,  or  Von  Hippel  1954,  page  64).  This  means 
the  wave  is  no  longer  linearly  polarized.  However,  in  the  case  of  normal  incidence, 
where  0;  is  equal  to  zero,  the  transmitted  and  reflected  waves  are  also  normal  to  the 
interface  with  no  polarization  changes. 

The  amplitudes  of  the  transmitted  and  reflected  waves  are  dependent  on  the 
incident  wave  amplitude,  the  material  properties  of  the  two  media,  and  the  angle  of 
incidence  (and  the  polarization  which  is  still  assumed  to  be  horizontal,  E,  or 
perpendicular  to  the  plane  of  the  incident  and  reflected  rays).  The  Fresnel  equations 
(Lorrain  and  Corson  1970,  page  508;  Kline  and  Kay  1965,  page  179)  describe  these 
relationships  and  the  formulas  are 

Er=  REj  and  E,=  TEj, 

where  R  and  T  are  the  reflection  and  transmission  coefQcients,  respectively,  and  are 
as  follows: 

COS0,.  -  u,  ^,COS0, 

R=-L2-1 — '  — J.  and 
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Note  that  for  complex  media,  these  coefficients  will  also  be  complex,  resulting  in 
amplitude  and  phase  changes  in  the  reflected  and  transmitted  waves. 
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3  Forward  Modeling  of  GPR 
Data 


Only  within  the  past  five  years  have  serious  efforts  been  made  to  model  full 
waveform  GPR  data.  These  efforts  have  concentrated  solely  on  the  forward 
modeling  process.  The  forward  modeling  routine  must  include  the  physical 
parameters  introduced  in  the  GPR  theory  section  (complex  relative  dielectric 
permittivity,  complex  relative  magnetic  permeability,  electrical  conductivity,  layer 
thickness,  and  time  relaxation  and  distribution  parameters)  and  model  the  effects 
these  parameters  have  on  the  transmitted  waveform.  The  processes  a  GPR  wave 
und^goes  as  it  passes  through  a  material  are  reflection  and  transmission  at 
interfaces  or  boundaries,  attenuation,  geometric  spreading,  dispersion,  and 
scattering.  The  GPR  forward  modeling  program  GPRMODv2  (Powers  and  Olhoeft 
1995)  employed  in  this  research  compensates  for  all  of  these  effects  except 
scattering.  The  program  makes  several  assumptions  to  simplify  the  modeling 
process;  zero  source-receiver  offset;  normal  incidence;  horizontal,  layered, 
homogeneous,  isotropic  media;  and  no  three  dimensional  electric  or  magnetic 
polarization,  near-field,  or  scattering  loss  effects.  The  assumption  of  zero  source- 
receiver  offset  and  normal  incidence  are  interrelated  and  will  be  discussed  later. 
Assuming  a  horizontally  layered  media  simplifies  the  wave  propagation  path  and  is 
a  reasonable  assumption  for  obtaining  an  initial  model  response  or  approximation  to 
field  data.  The  absence  of  electric  or  magnetic  polarization  is  a  direct  result  of 
assuming  normal  incidence  and  a  horizontal  interface.  In  reality,  the  electric  field 
probably  does  not  maintain  its  linearity  at  an  interface,  but  becomes  polarized  in 
some  elliptical  form.  However,  what  is  measured  at  the  receiver  is  only  the 
component  of  the  field  that  is  in  the  same  direction  as  the  transmitted  field. 

Reflection  and  transmission  of  the  GPR  signal  occur  at  boundaries  of  contrasting 
dielectric  properties.  With  each  boundary  encountered,  the  amplitude  of  the  original 
signal  is  attenuated  and  only  a  small  fi-action  of  the  originally  transmitted  signal  is 
measured  by  the  receiver,  hi  addition  to  the  primary  reflection  which  occurs  at  an 
interface  and  is  detected  by  the  receiver,  reflection  multiples  are  also  generated  at 
interfaces.  These  multiples  are  reflections  off  more  than  one  interface  as  the 
transmitted  signal  travels  through  the  subsurface.  It  is  important  to  differentiate  the 
primary  reflections  fi’om  multiple  reflections  since  multiples  do  not  represent  true 
interfaces  at  depth  (if  a  single  reflection  traveltime  is  assumed). 
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The  conductivity  of  a  material  is  a  major  factor  in  determining  the  degree  of 
signal  attenuation  The  higher  the  conductivity  value,  the  greater  the  amount  of 
signal  energy  that  is  absorbed  by  the  medium,  thus  reducing  the  amplitude  of  the 
signal  measured  at  the  receiver.  Soils  having  a  high  clay  and  moisture  content  can 
absorb  so  much  of  the  transmitted  signal  that  the  GPR  will  not  be  able  to  "see" 
below  that  layer.  If  identification  of  the  clay  layer  is  the  objective,  then  loss  of 
signal  at  depth  would  be  a  good  marker  for  identifying  the  clay.  However,  if  the 
layer  of  interest  is  below  the  clay  layo*,  then  another  technique  would  be  required  to 
locate  it 

There  are  two  types  of  geometric  spreading.  The  first  is  caused  by  the  antenna 
radiation  pattern  and  the  other  is  a  result  of  spherical  spreading  of  the  initial  energy 
ouqiut.  The  radiation  pattern  of  a  GPR  antenna  varies  the  amount  of  energy  output 
as  a  fimction  of  angle  from  the  normal  to  the  antenna.  In  addition,  the  radiation 
pattmi  varies  with  the  electromagnetic  properties  of  the  near  surface  material.  The 
coupling  effects  between  the  antenna  and  ground  surface  are  generally  ignored, 
however  the  forward  modeling  routine  used  in  this  research  does  provide  an  option 
for  fitting  the  model  response  to  the  field  data  by  compressing  or  stretching  the 
model  waveform.  The  effect  of  the  incident  angle  of  energy  output  on  the  reflection 
and  transmission  coefficients  was  investigated  using  the  Fresnel  equations  (chapter 
2).  Assuming  a  magnetic  permeability  of  free  space  and  positive,  real  propagation 
constants  that  increase  with  depth,  plots  of  the  reflection  coefficient  and 
transmission  coefficient  as  a  fimction  of  incident  angle  6  were  constructed  (Figure 
4).  It  can  be  seen  that  the  curves  are  relatively  flat  for  incident  angles  6<20°.  The 
effect  of  geometrical  spreading  caused  by  the  angle  of  incidence  of  the  energy  output 
is  related  to  the  separation  distance  between  the  transmitter  (Tx)  and  receiver  (Rx) 
antenna.  A  general  rule  of  thumb  given  by  Arman  and  Cosway  (1992)  states  that 
the  separation  between  the  transmitter  and  receiver  should  be  approximately  equal 
to  20%  of  the  target  depth  (Figure  5).  Based  on  this  assumption,  the  ratio  between 
one  half  the  Tx-Rx  distance  (s=2x)  and  the  target  depth  y  should  be  x/y=0. 1 .  This 
corresponds  to  an  incident  angle  of  0=:6°,  which  is  much  less  than  the  20°  limit. 
Thus,  if  the  general  rule  of  Arman  and  Cosway  (1992)  is  adhered  to,  normal 
incidence  can  be  assumed.  For  very  shallow  investigations,  high  frequency  anteimae 
are  required  which  have  very  small  source-receiver  offset  and  thus  virtually  normal 
incidence. 

The  second  type  of  geometric  spreading  occurs  when  the  transmitted  energy 
propagates  into  the  earth  and  expands  over  an  ever  increasing  spherical  surface. 

The  spreading  wavefront  is  only  truly  spherical  when  the  velocity  profile  of  the 
material  is  isotropic  and  homogeneous.  The  forward  modeling  program  corrects  for 
geometric  spreading  using  an  algorithm  based  on  May  and  Hron  (1978)  and  Balanis 
(1989). 

Dispersion  occurs  when  the  velocity  of  a  material  varies  with  frequency.  Both 
dielectric,  magnetic,  and  conductive  losses  can  cause  dispersion.  It  is  often  assumed 
that  the  velocity  of  a  layer  is  constant,  although  in  reality  the  velocity  is  generally 
frequency  dependent  over  at  least  some  frequency  range.  Mathematically, 
dispersion  is  represented  by  the  imaginary  components  of  the  dielectric  permittivity 
and  magnetic  permeability.  As  discussed  earlier,  the  complex  parameters  are 
introduced  via  the  Cole-Cole  model.  Dispersion  is  simulated  in  the  forward 
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modeling  program  by  computing  a  separate  time  delay  for  each  frequency  for  the 
traveltime  to  a  particular  reflection  event  (Powers  and  Olhoeft  1995). 

Scattering  occurs  on  a  small  scale  (scale  of  the  wavelength)  and  can  be  classified 
as  either  volume  scattering  or  surface  scattering.  Volume  scattering  is  caused  by 
objects  dispersed  in  the  media,  such  as  gravel  or  rocks,  which  scatter  the  signal  in 
multiple  directions  when  the  electromagnetic  wave  encounters  it.  Surface  scattering 
occurs  vdien  the  traveling  wave  hits  a  rough  interface  which  causes  the  wave  to 
spread  in  various  directions.  The  primary  effect  of  scattering  is  the  reduction  of 
signal  strength  (Olhoeft;  1984).  Neither  type  of  scattering  is  considered  in  the 
forward  modeling  program. 


Forward  Modeling  Example 

The  forward  modeling  process  begins  by  choosing  the  number  of  layers  in  the 
model  and  specifying  the  parameter  values  for  each  layer.  The  parameters  include 
layer  thickness,  electrical  conductivity,  relative  dielectric  permittivity,  relative 
magnetic  permeabilify,  and  the  Cole-Cole  relaxation  and  distribution  parameters  for 
the  permittivity  and  permeability.  Generally,  both  the  inexperienced  and 
experienced  modeler  do  not  have  an  intuitive  feeling  for  the  parameters  required  for 
modeling  GPR  data.  The  low  and  high  frequency  limits  for  both  the  real  relative 
dielectric  permittivity  and  real  relative  magnetic  permeability  are  required  as  input. 
The  default  values  for  the  magnetic  permeability  are  generally  suitable  for  modeling 
most  GPR  data;  only  in  cases  where  the  soil  has  a  significant  iron  content  will  these 
values  need  adjusting.  Values  for  the  relative  dielectric  permittivity  and  electrical 
conductivity  can  be  obtained  from  publications  of  laboratory  measurements  on 
similar  materials.  For  the  Cole-Cole  parameters,  the  default  values  are  good  initial 
values.  The  distribution  parameter  a  governs  the  slope  of  the  relaxation  curve 
whereas  the  relaxation  parameter  x  shifts  the  curve  laterally.  These  parameters  can 
be  determined  from  laboratory  measurements  or,  for  sand/clay  mixtures, 
GPRMODv2  includes  an  option  for  estimating  the  Cole-Cole  parameters  based  on 
the  sand/clay/water  ratio  (Olhoeft  1986). 

An  example  of  the  forward  modeling  output  is  shown  in  Figure  6.  The  model  is 
for  a  two  layer  rigid  pavement  structure  consisting  of  concrete  and  gravel  base 
overlaying  a  sandy  soil.  The  primary  reflections  A  and  B  in  the  output  trace  are 
caused  by  the  concrete/gravel  interface  at  0.2m  depth  and  the  gravel/sand  interface 
at  0.6m  depth,  respectively.  The  reflection  at  about  10ns  (a)  is  a  first  multiple  (i.e. 
wave  that  traveled  same  path  as  primary  reflection  A)  of  the  concrete/gravel 
interface,  vriiereas  reflection  b  is  the  first  multiple  of  the  gravel/sand  interface. 
Reflections  c  and  d  are  multiples  that  traveled  a  more  complex  path.  The  thickness 
of  layer  1  (concrete)  dictates  how  far  in  time  the  primary  reflection  A  is  shifted  from 
zero;  as  thickness  increases,  two-way  traveltime  increases.  The  separation  between 
primary  reflectors  A  and  B  is  controlled  by  the  thickness  of  the  gravel  base  (layer 
2).  As  this  thickness  decreases,  it  becomes  more  difficult  to  distinguish  the  location 
in  time  of  the  concrete/gravel  interface  (B)  because  the  onset  of  the  reflection  event 
gets  mixed  in  with  the  arrivals  of  the  primary  reflection  A  and  multiple  reflection  a. 
The  ability  to  differentiate  reflection  events  is  a  major  concern  in  the  inversion  of 
GPR  data. 
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Forward  modeling  is  a  usefiil  aid  in  determining  if  it  is  feasible  to  use  GPR  at  a 
particular  site,  what  frequency  antenna  is  required  to  detect  layers  of  varying 
thickness,  and  what  type  of  response  should  be  expected  for  various  combinations 
of  conductivity,  dielectric  permittivity,  and  magnetic  permeability.  Although 
forward  modeling  is  a  valuable  learning  aid,  one  is  often  in  the  situation  where  the 
GPR  data  has  been  collected  and  it  is  necessary  to  determine  the  layer  structure  and 
dielectric  properties  based  solely  on  the  radar  data;  this  is  when  the  inversion  of 
GPR  data  is  required.  The  inverse  modeling  of  GPR  data  is  discussed  in  Chapter  5. 
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4  Full  Waveform  GPR  Forward 
Modeling  Investigations 


Modeling  Problem  Selection  and  Details 

Problem  selection  and  justification 

The  concepts  of  EM  wave  propagation  theory  and  fiill  waveform  forward 
modeling  are  discussed  in  the  preceding  chapters.  In  this  chapter  a  common 
geotechnical  problem  is  select^  for  a  detailed  modeling  study  to  illustrate  the  use  of 
full  waveform  GPR  modeling  and  to  investigate  key  concepts  of  detectability  and 
resolution.  Included  in  the  modeling  study  are  effects  of  reflection  from  interfaces 
between  different  materials,  changes  in  propagation  path  length,  transmitter 
frequency,  and  material  properties,  multiple  reflection  events,  and  attenuation 
caused  by  finite  conductivity.  The  effects  of  dispersion  are  not  considered  in  the 
present  study  and  are  likely  not  significant  in  most  materials  over  the  frequency 
range  100  MHz  to  1  GHz  (Olhoeft  and  Capron  1994). 

Detection  of  air-  and  water-filled  cavities  in  concrete  structures  or  behind 
concrete  layers  is  the  problem  selected  for  study.  This  is  a  frequently  occurring 
problem  that  can  threaten  the  integrity  and  safe  operation  of  critical  structures  and 
an  application  for  which  GPR  is  commonly  considered.  Cavities  can  occur  in 
concrete  used  as  structural  components  (pillars,  walls  and  monoliths)  for  a  variety  of 
reasons,  e.g.,  as  an  artifact  of  the  construction  process,  a  result  of  concrete 
degradation,  and  as  a  result  of  delamination  or  debonding  from  other  structural 
components.  Concrete  is  also  used  in  the  form  of  layers  or  slabs  for  building  floors 
and  walls,  flow  channels,  soil  and  rock  slopes,  and  tunnel  liners.  Cavities  also  occur 
behind  concrete  layers  or  slabs  due  to  a  variety  of  processes,  such  as  differential 
settlement  or  piping.  Application  of  GPR  to  detection  of  cavities  within  or  behind 
concrete  is  always  a  tradeoff  between  the  requirement  for  high  frequencies  for 
resolution  of  cavity  thickness  or  of  thin  concrete  layers  and  lower  frequencies  for 
achieving  the  required  depth  of  investigation. 


Model  problem  sets 

Three  transmitter  center  frequencies,  300, 500,  and  900  MHz,  are  selected  for 
the  model  calculations,  and  are  representative  of  commercially  available  GPR 
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systems  and  typical  of  frequencies  which  might  be  selected  for  investigation  of 
concrete  structures.  The  Ricker  waveform  is  used  for  all  three  frequencies;  thus  the 
shape  of  the  pulse  is  the  same  for  all  three  frequencies,  with  only  the  duration 
varying  to  produce  the  variation  in  center  frequency.  One  major  GPR  manufacturer 
tnnintning  Aat  the  Ricker  waveform  is  an  extremely  close  approximation  to  the 
actual  waveform  for  their  antennae.  Duke  (1990)  investigated  the  anteimae  of  a 
second  major  GPR  manufacturer  and  produced  calibrated  transmittCT  waveforms; 
these  actual  calibrated  waveforms  are  compared  to  the  ideal  Ricker  waveforms  in 
Figure  7  for  the  same  designated,  nominal  center  frequencies.  Both  the  actual  and 
the  Ricker  waveforms  are  scaled  to  the  same  arbitrary  maximum  amplitude.  The 
actual  300  MHz  transmitter  waveform  and  spectra  are  very  similar  to  the 
corresponding  Ricker  waveform  and  spectra.  The  actual  500  and  900  MHz 
transmitter  waveforms  are  different  from  the  Ricker  waveforms  in  various  aspects, 
due  to  transmitter  ringing  and  noise. 

Two  problem  sets  are  considered  for  the  modeling  investigation.  Problem  Set  1 
consists  of  models  of  concrete  over  air  or  water.  For  the  concrete  over  air  cases,  the 
thickness  of  concrete  is  systematically  varied  from  0. 1  to  4  m,  while  for  the  concrete 
over  water  cases  a  fixed  concrete  thickness  is  selected.  Problem  Set  2  consists  of 
models  of  a  concrete  layer  over  a  water  layer  over  “sand”.  Cases  in  Problem  Set  2 
have  a  fixed  concrete  tUckness,  while  the  thickness  of  the  water  layer  is 
systematically  varied.  Also,  for  some  of  the  cases  in  Problem  Set  2,  GPR  forward 
models  are  conducted  with  and  without  multiple  events. 

Concrete  is  not  a  simple  material  to  characterize  generically.  The  EM  properties 
for  concrete  used  for  this  investigation  were  deduced  from  Buyukozturk  and  Rhim 
(1993)  for  three  concrete  water  contents  (w.c.),  9.5%,  6.7%,  and  0%,  as  shown  in 
Table  1,  which  also  includes  the  properties  for  air,  water  and  sand.  A  small 
frequency  dependence  of  the  relative  dielectric  permittivity  is  included  for  concrete. 
For  purposes  of  this  study,  this  representative  set  of  EM  properties  is  selected  based 
on  published  information,  without  consideration  of  the  relevance  of  specific  concrete 
mix  design,  i.e.,  type  of  cement,  sand/gravel/cement/water  mix,  type  and  size 
distribution  of  aggregate,  etc.  The  concrete  is  considered  massive  and  unreinforced. 
The  concrete  specimens  used  for  the  measurements  by  Buyukozturk  and  Rhim  were 
cast  with  a  water/cement/sand  mix  ratio  of  1 :2:4  by  weight,  with  no  coarse 
aggregate.  AportlandcementofTypelwasused.  Specimen  age  was  4  weeks  at 
the  time  of  measurements.  The  effect  reinforced  concrete  has  on  the  penetration  of  a 
transmitted  GPR  signal  was  not  studied  in  this  investigation. 
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Table  1 

EM  Properties  of  Model  Materials 

Model  Material 

Relative  Dielectric 
Permittivity  (Real) 

Conductivity 

mS/m 

300 

500 

MHz 

900 

Concrete  w.c.  =  9.5% 

8.4 

8.3 

8.1 

40 

w.c.  =  6.7% 

5.2 

5.1 

5.0 

30 

W.C.  =  0.0% 

3.5 

3.55 

3.6 

10 

Air 

1 

1 

1 

0 

Water 

80 

80 

80 

20 

Sand 

25 

25 

25 

5 

Conceptual  modeling  considerations 

A  conceptual  understanding  of  the  GPR  modeling  process  requires  consideration 
of  EM  wave  propagation  velocities,  two-way  reflection  times  from  interfaces  in  the 
model,  reflection  and  transmission  coefficients,  effects  of  geometric  spreading  and 
attenuation  on  wave  amplitudes,  and  multiple  reflection  events.  An  imderstanding 
of  the  modeling  process  also  is  extremely  helpful  if  not  necessary  for  formulating  or 
deducing  a  subsurface  model  that  explains  or  is  consistent  with  the  features 
observed  in  field  GPR  data.  To  illustrate  the  process  of  conceptual  modeling,  the 
model  used  in  Problem  Set  2,  concrete  over  a  water-filled  cavity  over  sand,  is 
utilized. 

Consider  the  specific  case  of  a  0.3  m  thick  concrete  layer  over  a  0.05  m  thick 
water  cavity  over  a  saturated  sand,  which  is  considered  infinitely  thick  for  purposes 
of  the  modeling  investigation.  The  EM  wave  propagation  velocity  (speed)  C;  and 
EM  wavelength  A.; ,  where  the  subscript  i  refers  to  the  i-th  material  in  the  model, 
and  the  amplitude  reflection  coefficients  R^j  and  transmission  coefficients  T^ , 
where  ij  refers  to  the  interface  between  the  i-th  and  j-th  materials  in  the  model,  are 
given  by  the  following  equations 

Ci  “  Cq/  Eji'  , 

Ry  =  (v/e;  -  )  /  (ycri'  +  ) , 

and 

Tij  =  276,7  (v/e;  +  ye,,'), 

\siiere  Cq  is  the  EM  wave  speed  in  a  vacuum  (3x10*  m/s),  e,'  is  the  real  part  of  the 
relative  dielectric  permittivity  of  the  i-th  model  material,  and /is  the  frequency.  The 


18 


Chapter  4  Full  Waveform  GPR  Forward  Modeling  Investigations 


above  expressions  for  the  reflection  and  transmission  coefficients  are 
simplifications  for  the  full  fi-equency  dependent  expressions  (see  Chapter  2)  used  in 
GPRMODv2.  Using  the  subscripts  c,  w,  and  s  for  concrete,  water  and  sand, 
respectively.  Table  2  gives  the  values  for  wave  speed,  wavelength,  reflection 
coeflBcient,  and  transmission  coefficient  for  the  model  materials.  The  velocities  for 
concrete  are  averages  for  the  three  fi’equencies  for  each  water  content  value. 


Table  2 

EM  Velocities,  Wavelengths,  Reflection  and  Transmission 
Coefficients  _ 


Model  Material 


Velocity, 
xK)*  nVs(m/ns) 


Wavelength,  m 
300  500  900 


Reflection 

Coefficient 


Transmission 

Coefficient 


Layer  1: 

Concrete 
w.c.=  9.5% 
w.c.=  6.7% 
w.c.=  0.0% 


1.04  (0.104) 
1.33  (0.133) 
1.59  (0.159) 


.35  .21  .12 
.44  .27  .15 
.53  .32  .18 


■  0.515 
■0.599 
-0.654 


0.485 

0.401 

0.346 


Layer  2: 
Water 


0.33  (0.033) 


.11  .07  .04 


Layer  3: 
Sand 


0.6  (0.06) 


.20  .12  .08 


With  the  wave  velocities  in  Table  2,  it  is  possible  to  determine  the  two-way 
travel  times  through  the  layers  of  the  model.  Figure  8  shows  the  results  of  the 
traveltime  calculations  for  selected  GPR  events  in  the  model,  where  the  terminology 
“event”  refers  to  the  arrival  of  energy  at  the  receiver  that  has  followed  a  distinct 
propagation  path  through  the  model  (or  the  subsurface  in  the  case  of  actual  GPR 
data).  Thus  the  GPR  records,  either  field  data  or  model  calculations,  consist  of 
much  more  than  just  the  primary  reflection  events  from  each  interface.  In  some 
cases  identification  of  the  multiple  events  is  helpful,  but  in  most  cases  the  multiples 
complicate  the  interpretation  process  by  interfering  (superimposing)  with  both  the 
primary  reflections  and  with  other  multiples.  It  is  also  possible  that  interfaces  will 
be  so  close  together  (i.e.,  a  thin  layer)  that  the  primary  reflections  will  interfere  with 
each  other  and  the  two  interfaces  will  not  be  resolved.  The  numerical  designation 
of  the  multiples  refers  to  the  interfaces  where  reflections  occur  along  the  multiple 
propagation  path,  where  the  top  or  air/layer  1  interface  is  assigned  ‘O’.  Thus  the 
designation  ‘101’  refers  to  a  event  which  reflects  from  interface  1,  interface  0,  and 
interface  1;  this  multiple  event  is  shown  as  the  third  event  from  the  top  in  Figure  8. 

Whether  or  not  events  interfere  with  each  other  or  a  layer  is  resolved  depends  on 
the  frequency  and  hence  the  duration  and  dominant  wavelength  of  the  propagating 
waveform.  TTie  dominant  wavelengths  of  the  nominal  Ricker  waveforms  are  shown 
in  Table  2.  Widess  (1973)  demonstrates  that  a  layer  is  theoretically  resolved  for 
thicknesses  greater  that  Vi,  where  k  is  the  dominant  wavelength  of  the  waveform 
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within  the  layer.  However,  practical  resolution  is  typically  taken  to  occur  at 
approximately  A/4  (Jol  1995),  and  “visual  resolution,”  i.e.,  the  thickness  for  which 
the  two  primary  reflections  are  visually  obvious,  occurs  at  approximately  A/2.  Thus 
in  the  absence  of  interfering  multiples,  visual  resolution  of  the  water  layer  (see 
Table  2)  should  occur  for  layer  thickness  of  approximately  0.06, 0.04,  and  0.02  m 
for  frequencies  of  300, 500,  and  900  MHz,  respectively.  Based  on  these 
considerations,  the  water  layer  in  Figure  8  (thickness  hw  =  0.05  m)  should  be 
visually  resolved  for  500  and  900  MHz  but  not  for  300  MHz. 

A  useful  way  to  visualize  a  GPR  record  is  as  the  convolution  of  the  propagating 
waveform  with  the  impulse  response  of  a  model  or  the  earth.  The  impulse  response 
is  constructed  simply  as  a  spike  at  the  two-way  reflection  time  for  each  event  with 
amplitude  and  polarity  equ^  to  the  reflection  coefficient  (see  Table  2);  this 
simplistic  impulse  response  ignores  transmission  and  propagation  losses.  In  the 
convolutional  model  response,  each  event  will  have  duration  equal  to  duration  of  the 
propagating  waveform.  Interference  will  then  occur  when  events  occur  closer  in 
time  than  the  duration  of  the  propagating  waveform.  From  Figure  7,  a  duration  for 
the  Ricker  waveforms  can  be  defined  as  the  visually  apparent  time  between  zero 
departure  (total  duration)  for  the  waveform  or  as  the  time  between  the  negative 
peaks  of  the  waveform.  The  approximate  values  for  these  definitions  of  duration 
are  shown  below; 


300  MHz  500  MHz  900  MHz 
Total  Duration,  ns  6.3  3.7  2.15 

Time  between  negative  peaks,  ns  2.5  1.5  0.85 

Comparing  these  values  for  duration  to  the  two  primary  reflection  two-way  times  in 
Figure  8  indicates  that  the  two  primary  events  should  be  well  resolved  for  900  MHz, 
interfere  slightly  for  500  MHz  (although  the  two  negative  peaks  for  each  event 
should  be  resolved),  and  should  significantly  interfere  for  300  MHz.  The  resolution 
assessment  based  on  waveform  duration  is  completely  consistent  with  the 
assessment  based  on  dominant  wavelength  in  the  preceding  paragraph. 

Results  of  GPRMODV2  calculations  for  the  concrete/water/sand  model  of 
Figure  8  with  concrete  w.c.  =  9.5%  and  900  MHz  Ricker  waveform  are  shown  in 
Figure  9;  the  top  plot  is  for  no  multiples,  and  the  bottom  plot  includes  multiple 
events.  The  two  primary  events  are  labeled  and  are  clearly  resolved  as  predicted 
previously  for  900  MHz  based  on  waveform  duration  and  dominant  w'avelength 
considerations.  Two  multiple  events  are  indicated  in  the  bottom  plot  The  first 
labeled  multiple  event  is  actually  the  superposition  of  two  multiple  events,  101  and 
212,  that  arrive  nearly  simultaneously  (see  Figure  8  under  w.c.  =  9.5%).  The 
multiple  event  201  is  attenuated  significantly,  and  all  other  multiple  events  are  not 
evident. 

As  an  illustration  of  the  effects  of  propagation  velocity  on  the  spacing  of  primary 
and  multiple  events.  Figure  10  shows  GPRMODv2  calculation  results  for  concrete 
w.c.  =  0.0%  and  900  MHz,  i.e.,  identical  to  the  model  for  Figure  9  except  for  water 
content  and  hence  propagation  velocity  (see  Table  2).  From  the  reflection  timesnn 
Figure  8  under  w.c.  =  0.0%  and  dominant  wavelength  and  duration  considerations,  it 
is  evident  that  the  primaries  should  still  be  resolved;  however,  multiple  event  101 
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should  now  interfere  with  the  2nd  primaiy,  while  multiple  events  101  and  212 
should  be  resolved.  All  of  these  predictions  are  evident  in  the  calculated  GPR 
model  response  and  the  events  are  identified  in  Figure  10.  Multiple  event  202  is 
also  evident  and  identified  for  this  case,  because  of  the  decreased  attoiuation  for  the 
lower  conductivity  associated  with  w.c.  =  0.0%  compared  to  the  higher  conductivity 
with  w.c.  =  9.5%. 


As  a  final  example  to  illustrate  forward  modeling  concepts,  the  GPRMODv2 
results  for  the  previous  model  with  w.c.  =  9.5%  and  using  a  300  MHz  propagating 
waveform  are  shown  in  Figure  11.  Considering  the  arrival  times  in  Figure  8  under 
w.c.  =  9.5%  and  dominant  wavelength  and  duration  for  300  MHz,  it  is  obvious  that 
not  only  should  the  two  primary  events  not  be  resolved  but  multiple  events  101  and 
212  should  interfere  with  the  primaries.  The  GPR  model  results  in  Figure  1 1 
illustrate  the  lack  of  resolution  of  the  two  primaries  and  comparison  of  the 
calculation  wdth  no  multiples  (top  plot)  with  the  calculation  with  multiples  (bottom 
plot)  illustrates  the  interference  of  multiple  events  101  and  212  with  the  primaries. 
None  of  the  other  multiples  are  evident,  illustrating  greater  attenuation  of  the  300 
MHz  waveform  relative  to  the  900  MHz  case  considered  earlier.  Although  the 
complex  waveforms  in  Figure  1 1  suggest  superposition  of  events,  arrival  times 
other  than  that  of  the  1st  primary  cannot  be  determined  with  confidence  by  visual 
inspection.  Application  of  data  processing  procedures,  such  as  crosscorrelation  and 
deconvolution,  to  detect  event  arrivals  that  are  too  closely  spaced  in  time  to  be 
visually  resolved  is  discussed  in  the  next  chapter. 


Scope  and  Results  of  the  Modeling  Study 

Scope  of  study 

Appendix  A  contains  tabulations  of  the  cases  included  in  the  Problem  Sets. 
Table  A 1  lists  the  27  cases  of  Problem  Set  1  for  concrete  over  air  for  which  GPR 
model  responses  were  computed  and  reproduced  for  this  report.  Table  A2  lists  the 
cases  of  Problem  Set  1  for  concrete  over  water.  The  GPR  model  response 
calculations  fi-om  GPRMODv2  for  Problem  Set  1  are  contained  in  Appendix  B. 
Table  A3  lists  the  cases  of  Problem  Set  2  for  models  of  concrete  over  water  over 
sand .  Cases  37-144  in  Table  A3  are  for  calculations  with  multiples,  while  cases 
145-252  are  for  calculations  with  no  multiples.  The  GPR  model  response 
calculations  for  Problem  Set  2  are  contained  in  Appendix  C. 


Selected  results  and  discussion:  Problem  set  1 

Results  of  the  GPRMODv2  model  calculations  for  Problem  Set  1  are  presented 
in  Appendix  B.  Figure  12  combines  Figures  B4  and  B1 1  for  a  side-by-side 
comparison  example.  The  three  plots  on  the  left  side  of  Figure  12  are  for  concrete 
over  air,  while  the  three  plots  on  the  right  are  for  concrete  over  water,  with  the  same 
fi'equencies  side-by-side  and  concrete  thickness  h  =  0.3  m  and  w.c.=9.5%  for  all' 
cases.  All  of  the  concrete/air  and  concrete/water  cases  show  the  primary  and  first 
multiple  event.  From  Table  2,  the  reflection  coefficient  for  concrete  over  water  is 
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-0.5 15,  while  the  reflection  coefiQcient  for  concrete  over  air  is  +0.485.  Thus  the  first 
expected  and  observed  difference  in  the  concrete/air  and  concrete/water  cases  are  the 
diffwftnt  polarities  of  the  reflection  events.  For  the  concrete  thickness  0.3  m,  the 
primaries  and  multiples  are  well  separated  for  all  three  frequencies.  The  concrete 
thickness  for  this  problem  can  be  determined  from  the  model  response  by  using  the 
arrival  time  of  the  primary  or  the  time  difference  between  primary  and  multiple 
multiplied  by  the  propagation  velocity.  The  approximate  two-way  time  diffCTence 
between  primaries  and  multiples  is  t  =  6  ns;  the  concrete  thickness  h  is  given  by 
h  =  (t/2)(c/ye,), 

h  =  0.5  X  (6  X  10  ®)  X  (3  X  10*)// 8.4  m  =  0.31  m, 
in  satisfactory  agreement  with  the  known  model  value. 

Another  more  subtle  feature  of  the  cases  in  Figure  12  is  the  effect  on  amplitudes 
of  the  reflection  process  and  total  signal  attenuation.  For  each  plot,  a  “Display 
Gain”  is  shown  (rounded  to  the  nearest  dB).  This  gain  is  a  constant  by  which  all 
values  in  the  model  calculation  output  values  are  multiplied  to  make  the  maximum 
amplitude  in  the  output  equal  to  80%  of  the  maximum  value  for  the  graphic  display 
(Powers  and  Olhoeft  1995).  Thus  the  Display  Gain  is  directly  relate!  to  the  total 
signal  loss  as  a  result  of  attenuation,  spreading,  and  reflection  and  transmission  at 
boundaries  for  all  cases.  In  fact,  the  total  signal  loss  will  always  be  greater  than  the 
display  gain  in  dB,  since  the  display  gain  will  be  determined  by  the  amplitude  of  the 
first  primary  event  rather  than  the  amplitude  of  the  initial  transmitter  waveform  (the 
first  primary  amplitude  will  always  be  less  than  the  initial  transmitter  waveform 
ampUtude).  Two  features  are  apparent  from  Figure  12:  (1)  signal  loss  increases  as 
frequency  increases  for  both  primaries  and  multiples  as  expected;  (2)  signal  losses 
are  greater  by  approximately  1  dB  for  the  concrete/air  cases  than  for  the 
concrete/water  cases.  The  second  feature  is  understood  by  comparing  the  reflection 
coefBcients  for  concrete/air  and  concrete/water.  Only  48.5%  of  the  incident 
amplitude  is  reflected  at  the  air  boundary,  while  5 1.5%  is  reflected  at  the  water 
boundary.  This  difference  in  reflected  amplitude  at  the  two  boundaries  amounts  to 
IdB  (  by  converting  the  reflection  coefficients  to  dB  and  comparing)  as  shown  by 
the  model  calculations. 

In  addition  to  the  cases  shown  in  Appendix  B,  for  concrete  thicknesses  of  0. 1, 

0.3,  and  0.75  m,  the  thickness  h  was  systematically  varied  from  0.01  m  to  4.0  m  to 
determine  signal  loss  as  a  function  of  thickness  for  the  GPR  model  calculations. 
Signal  loss  in  dB  versus  concrete  thickness  for  concrete  w.c.  =  9.5%  (top)  and  0.0% 
(bottom)  is  plotted  in  Figure  13,  where  the  “signal  loss”  is  just  set  numerically 
equal  to  the  display  gain.  The  signal  loss  plots  illustrate  several  important  concepts 
regarding  the  interplay  of  electromagnetic  parameters  in  determining  total  signal 
loss  and  achievable  penetration  depth  in  concrete.  The  plots  in  Figure  13  have  the 
appearance  of  two  linear  regions  joined  by  a  transition  region.  The  initial  steeply 
sloping  region  of  the  plots  corresponds  to  the  region  where  geometric  spreading 
dominates  the  total  signal  loss,  while  the  second  lesser  sloping  region  corresponds 
to  domination  of  the  total  signal  loss  by  actual  intrinsic  propagation  loss.  Wiy  the 
transition  region  is  so  abrupt,  particularly  for  the  w.c.  =  9.5%  cases,  is  not  clear.' 

The  initial  linear  regions  for  the  w.c.  =  9.5%  cases  and  for  the  w.c.  =  0.0%  cases  are 
practically  coincident,  while  the  slope  for  the  w.c.  =  9.5%  cases  is  steeper  than  for 
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the  w.c.  =  0.0%  cases.  The  slopes  within  the  second  linear  region  for  the  three 
frequencies  are  approximately  equal  for  each  w.c.;  although  the  curves  for  each 
frequency  are  separated  by  approximately  5  dB,  i.e.,  signal  loss  for  500  MHz  is  5 
dB  greater  than  300  MHz  and  signal  loss  for  900  MHz  is  5  dB  greater  than  500 
MHz.  The  second  linear  region  begins  at  approximately  0.75  m  for  w.c.  =  9.5%  and 
1.75  m  for  w.c.  =  0.0%.  Slopes  within  the  second  linear  region  approximately 
equal  1.4  dB/m  for  concrete  w.c.  =  9.5%  and  approximately  equal  1  dB/m  for 
concrete  w.c.  =  0.0%,  consistent  with  observed  values  for  similar  materials  reported 
by  others  (Annan  1992). 

A  somev^at  initially  surprising  result  is  that  the  plots  predict  greater  signal  loss 
at  larger  propagation  distances  (>  2  m)  for  the  lower  conductivity  concrete.  The 
explanation  for  this  somewhat  counterintuitive  result  is  the  role  played  by  dielectric 
permittivity  in  signal  loss  and  the  geometric  spreadmg  formulation  (see  Chapter  3 
and  Powers  and  Olhoeft  1995).  The  factor  of  two  decrease  in  dielectric  permittivity 
plays  a  largo"  role  in  determining  penetration  than  the  factor  of  four  decrease  in 
conductivity;  this  relative  effect  may  not  be  the  same  for  other  ranges  of  permittivity 
and  conductivity. 

A  typical  value  for  total  loss  that  can  occur  from  all  propagation  effects  and  still 
result  in  detection  of  the  signal  at  the  Rx  is  100  dB  (Annan  1992).  Assuming  that 
the  signal  loss  plotted  in  Figure  13  is  the  actual  total  propagation  loss,  this  will  give 
a  conservative  estimate  of  achievable  maximum  concrete  thickness  penetration  for 
each  frequency  at  the  thickness  where  each  curve  crosses  the  100  dB  loss  value.  It 
is  assumed  that  the  average  of  the  maximum  thicknesses  for  w.c.  =  9.5%  and  0.0% 
is  the  tnavimniTi  thickness  for  an  “average  concrete”;  these  values  are  tabulated 
below: 


300  Mhz  500  Mhz  900  MHz 

Maximum  Penetration  Thickness,  m  8  6.8  3.7 

(For  ‘average  concrete”  at  100  dB 
total  signal  loss) 

These  values  of  100  dB  penetration  thickness  are  consistent  with  values  of  100  dB 
penetration  range  calculations  reported  elsewhere  for  similar  conductivities  and 
permittivities  (Annan  1992). 


Selected  results  and  discussion:  Problem  Set  2 

Problem  Set  2  is  a  systematic  investigation  of  the  primary  and  multiple  events 
associated  with  the  concrete/water/sand  model  discussed  in  the  Conceptual 
modeling  considerations”  section.  For  this  investigation,  the  water  layer  thickness  is 
varied  from  0.01  m  to  0.4  m  for  three  concrete  water  contents  and  three  frequencies. 
The  concrete  thickness  is  fixed  at  0.3  m  for  all  cases  in  Problem  Set  2.  The 
complete  GPRMODv2  calculation  results  of  the  investigation  are  presented  in 
Appendix  C  which  serves  as  a  reference  catalog  for  this  type  problem.  For 
problems  of  this  type,  the  key  factors  to  be  determined  are  concrete  thickness,  cavity 
thickness,  and  whether  the  cavity  is  water-  or  air-filled.  Determination  of  the 
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concrete  thickness  was  discussed  in  the  preceding  section.  Determination  of  the 
cavity  (water  layer)  thickness  requires  the  arrival  time  of  the  second  primary  event. 
Thus  the  two  primary  events  associated  with  the  cavity  must  be  resolved  either 
visually  or  via  data  processing.  Identification  of  the  two  primary  events  allows  the 
determination  of  the  concrete  and  cavity  thicknesses.  The  propagation  velocities 
within  the  concrete  and  cavity  are  also  determined,  and  in  turn  the  real  part  of  the 
dielectric  permittivity  of  the  concrete  and  cavity-filling  material,  which  will  identify 
the  cavity-fill  material  (air  or  water).  The  preceding  discussion  is  a  description  of  a 
manual  inversion  of  GPR  data.  In  the  actual  case  of  field  GPR  data,  all  the  multiple 
events  and  noise  will  be  present  to  complicate  the  identification  of  the  second 
primary  and  determination  of  its  arrival  time.  The  value  of  the  full  waveform 
forward  modeling  capability  is  for  iterative  manual  inversion  or  as  the  core  of  an 
automated  inversion  capability,  since  all  primary  and  multiple  events  (as  well  as 
attenuation,  spreading  and  dispersion  effects)  are  included. 

The  cases  for  water  layer  thickness  hw  =  0.05  m  are  discussed  in  detail  in  the 
“Conceptual  modeling  considerations”  section.  For  this  water  layer  thickness,  the 
two  primary  reflections  are  just  resolved  for  900  MHz  but  not  resolved  for  300 
MHz  (see  Figures  10  and  1 1).  Examination  of  the  results  in  Appendix  C  reveals  the 
following  approximate  water  layer  thicknesses  required  for  visual  resolution: 

300  MHz  500  MHz  900  MHz 

Minimum  Resolvable  Layer  Thickness,  m  0.1  0.05  0.03 

Although  there  is  a  slight  dependence  on  concrete  properties  affecting  the 
resolution,  the  above  visually  resolvable  thicknesses  are  good  rules-of-thumb. 

Considering  the  minimum  water  layer  thickness  required  for  visual  resolution,  as 
in  the  preceding  paragraph,  results  in  output  fi'om  GPRMODv2  in  which  the 
primary  and  multiple  events  are  closely  spaced  in  time  and  superimposed  to  some 
extent.  It  is  useful  and  instructional  to  consider  a  water  layer  thickness  for  which  all 
the  events  are  separated.  Figure  14  reproduces  Figure  C54,  Case  252,  for  water 
layer  thickness  hw  =  0.4  m.  An  additional  display  gain  of  20  dB  is  applied  to  Figure 
14  to  enhance  recognition  of  the  events.  Eight  events  are  identified  in  Figure  14. 
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5  GPR  Inverse  Modeling 


The  numerical  inversion  of  any  geophysical  data  is  basically  an  iterative  process 
where  a  forward  model  is  fit  to  the  data  and  the  parameters  adjusted  imtil  an 
acceptable  fit  between  the  model  and  the  data  is  achieved.  The  greater  the  number 
of  parameters  involved  in  the  inversion,  the  more  complex  the  inversion  process. 
Since  some  GPR  parameters  are  interrelated,  a  unique  inversion  solution  cannot  be 
obtained.  Because  of  the  large  number  of  parameters  involved  and  the 
nonuniqueness  of  the  solution,  the  usual  techniques  for  inverting  geophysical  data  to 
obtain  a  global  solution  of  all  parameters  simultaneously  via  an  iterative  least 
squares  inversion  procedure  (Simms  and  Morgan  1992)  is  not  practical  for  full 
waveform  GPR  data.  However,  experience  and  knowledge  of  typical  values  of 
parameters  for  a  given  material  can  aid  in  simpUfying  the  GPR  inversion 
requirements  and  in  determining  if  the  inversion  result  is  reasonable. 


Parameter  Considerations 

Based  on  the  forward  modeling  program  GPRMODv2  (Powers  and  Olhoeft 
1995)  discussed  previously,  a  minimum  of  thirteen  parameters  are  involved  in  the 
modeling  of  GPR  data.  Thus,  a  minimum  of  thirteen  parameters  would  be  required 
as  input  to  the  inversion.  The  parameters  include  three  that  pertain  to  the  overall 
model  (number  of  layers,  coupling  ratio,  and  zero  time  offset)  and  ten  that  are 
required  for  each  layer  (conductivity,  thickness,  and  eight  Cole-Cole  parameters). 
The  number  of  initial  variable  inversion  parameters  can  be  reduced  to  four  if  some 
simplifying  assumptions  are  made.  The  parameter  zero  time  of&et  is  obtained  from 
header  information  when  the  data  were  collected  in  the  field;  it  is  dependent  on  the 
GPR  system  used  to  acquire  the  data.  A  second  parameter,  coupling  ratio,  is 
determined  from  the  expression  This  relationship  was  obtained  by 

assuming  the  antenna  is  designed  for  a  material  having  a  dielectric  permittivity  of 
e,=4  so  the  velocity  =/doig„^„tenn.  =  c/Ve,  =  c/2.  The  interaction  of  the 
transmitted  pulse  with  the  ground  surface  tends  to  reduce  the  peak  operating 
frequency  of  the  antenna  thus  =  (/Lign^ant*imJ/Cr=  c/ZCr.  A  permittivity  of 
four  is  arbitrary  but  the  expression  for  the  coupling  ratio  provides  an  adequate 
initial  approximation.  The  two  parameters  dealing  with  the  magnetic  permeability 
can  assumed  initially  to  be  constant,  their  values  governed  by  whether  the  material 
is  or  is  not  magnetic.  However,  if  a  material  is  thought  to  have  a  significant  iron 
ccmtent  then  these  parameters  should  be  allowed  to  vary  throughout  the  inversion. 
The  four  Cole-Cole  relaxation  and  distribution  parameters,  two  related  to  the 
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dielectric  permittivity  and  two  to  the  magnetic  permeability,  can  also  be  initially 
taken  as  constant.  If  the  material  is  dispersive  or  wet,  then  the  default  values  will 
need  to  be  adjusted.  The  five  remaining  parameters,  number  of  layers,  layer 
thickness,  electrical  conductivity,  and  hi^i/low  relative  dielectric  permittivities,  are 
the  major  parameters  involved  in  modeling  GPR  data.  For  the  inversion  though,  the 
number  of  layers  and  layer  conductivity  (if  measured  in  the  field  or  laboratory)  can 
be  fixed,  \duch  leaves  three  parameters  per  layer  that  will  vary  throughout  the 
inversion  process,  layer  thickness  and  high/low  relative  dielectric  permittivities,  plus 
the  coupling  ratio. 

The  first  step  in  the  inversion  process  is  to  determine  how  many  layers  are 
represented  in  the  GPR  trace.  Two  methods  were  investigated  to  do  this, 
crosscorrelation  and  deconvolution.  Crosscorrelation  compares  two  sequences  to 
determine  how  similar  one  sequence  is  to  the  other.  An  example  of  using 
crosscorrelation  to  identify  interface  boundaries  is  shown  in  Figure  15.  In  Figure 
15a,  a  900  MHz  Ricker  waveform  was  crosscorrelated  with  the  forward  model 
resjjonse  from  a  three  layer  model  consisting  of  concrete,  water,  and  sand  (no 
multiples).  A  Ricker  waveform  is  representative  of  the  output  signal  for  a 
commercially  available  GPR  antenna  system.  The  top  plot  is  the  Ricker  waveform, 
the  center  plot  is  the  model  response,  and  the  bottom  plot  is  the  result  of  the 
crosscorrelation.  Figure  15b  represents  the  same  set  of  plots  except  the  input 
wavelet  (top  plot)  is  that  for  a  900  MHz  antenna  of  another  GPR  manufacturer.  For 
both  input  waveforms,  the  crosscorrelation  results  in  a  replica  of  the  input  wavelet 
at  each  layer  interface.  The  apparent  discrepancy  in  the  location  of  the  interfaces 
relative  to  time  (location  of  peaks  on  the  model  response  and  crosscorrelation 
curves)  represents  the  time  difference  between  the  peak  on  the  input  wavelet  and 
zero  time.  The  crosscorrelation  ignores  any  time  offset  of  the  input  signal,  thus  time 
t'  of  peak  B'  on  the  crosscorrelation  curve  is  equal  to  time  t  between  peak  A  on  the 
input  wavelet  and  peak  B  on  the  model  response  curve. 

Deconvolution  was  also  tested  as  a  method  to  identify  layer  interfaces. 
Deconvolution  is  a  process  designed  to  restore  a  waveshape  to  the  form  it  had 
before  it  underwent  a  linear  filtering  action  (Sheriff  1984).  The  same  input  and 
model  response  curves  used  above  were  deconvolved  to  obtain  the  location  of  the 
layer  boundaries  (Figure  16).  Again,  the  input  signal  in  Figure  16a  is  a  Ricker 
waveform  and  Figure  16b  is  the  other  calibrated  antenna  waveform.  The 
deconvolution  is  also  relative  to  the  input  wavelet  and  ignores  any  initial  time  offset. 
Note  that  the  deconvolution  output  is  similar  to  that  of  an  impulse  rather  than  a 
replica  of  the  input  wavelet,  as  was  the  case  in  the  crosscorrelation. 

When  the  Ricker  waveform  is  used  as  the  input  signal,  both  the  crosscorrelation 
and  deconvolution  provide  a  relatively  sharp  response.  However,  as  the  input  signal 
becomes  more  complex,  as  with  the  other  input  waveform,  deconvolution  provides  a 
cleaner  and  sharper  identification  of  the  layer  interface.  The  simplification  of  the 
output  to  an  impulse-like  response  is  desirable  in  order  to  obtain  a  clearer  definition 
of  the  layer  boundary  and,  in  cases  when  thin  layers  are  present,  for  separation 
(resolution)  of  the  layers. 
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A  rou^  estimate  of  the  interface  depth  can  be  determined  based  on  an  estimate 
of  the  material  velocity  and  the  traveltime.  This  should  provide  an  adequate  initial 
value  fw  the  inversion  routine. 

The  electrical  conductivity  can  be  obtained  from  a  reference,  which  typically 
would  rqx)rt  a  range  of  values  for  a  given  material,  laboratory  measurements,  or 
data  coUected  in  the  field.  However,  it  is  preferable  to  measure  the  conductivity  of 
the  matoial  since  conductivity  values  can  vary  over  a  wide  range  and  are  dependent 
on  site  specific  conditions. 

Values  for  the  relative  dielectric  permittivity  are  obtained  from  laboratory 
measurranents  performed  on  field  samples  or,  more  commonly,  from  publications  of 
laboratory  measurements  on  similar  materials.  The  permittivity  of  the  near  surface 
soil  layCT  can  also  be  obtained  using  time  domain  reflectometry  (TDR).  This  field 
measuronent  would  correspond  to  the  high  frequency  dielectric  permittivity  in  the 
forward  modeling  program.  The  real  dielectric  permittivity  tends  to  range  between 
4  and  25  for  typicd  soils  and  pavements.  In  the  higher  frequency  range,/>100 
MHz,  the  high  frequency  dielectric  permittivity  dominates  and  the  low  frequency 
dielectric  permittivity  can  be  set  equal  to  However,  if  dispersion  is  a  factor 
within  the  operating  frequency  range  of  the  antenna,  then  should  be  set  greater 

than  e^.. 

A  list  of  the  thirteen  inversion  parameters  and  their  general  range  of  values  is 
given  below.  The  values  for  the  relative  dielectric  permittivity  are  based  on  one 
author's  experience  whereas  the  values  for  relative  magnetic  permeability  and  the 
Cole-Cole  parameters  were  obtained  from  forward  modeling  exercises  using 
GPRMODv2  (Powers  and  Olhoeft  1995).  Many  of  these  parameter  ranges  are 
based  on  laboratory  measurements  of  field  samples  (Olhoeft  1986,  Olhoeft  and 
Capron  1994).  Note  that  the  dielectric  permittivity  and  magnetic  i>ermeability 
parameter  values  have  -1  flags  in  the  program  to  indicate  a  metallic  or  magnetic 
object. 


PARAMETER  RANGE 


Number  of  Layers 
Layer  Thickness,  h 
Conductivity,  o 


•Ce 

Pri 


M 

Cr 

Offset 


variable 
variable 
0-50  mS/m 

4-25,  dry  sand  4-10,  metallic  object  -1 
dry  1.0,  wet  10.0,  iron  in  soil:  dry  0.8,  wet  10.0  ns 
dry  1.0,  wet  0.8,  iron  in  soil:  dry  0.7,  wet  0.8 
iron-free  soil  1.0,  iron  in  soil  1-16,  magnetic  object  -1 
iron-free  soil  1.0,  iron  in  soil  1-2,  magnetic  object  -1 
iron-free  soil  1.0,  iron  in  soil  >0.3  ns 
iron-free  soil  1.0,  iron  in  soil  0.75 
1-2 

dependent  on  antenna  system 
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After  the  initial  values  of  the  parameters  have  been  chosen,  the  inversion 
proceeds  by  calculating  a  model  GPR  trace  and  comparing  it  to  the  field  trace. 
Usually  some  measure  of  fit  is  used  to  determine  if  the  model  trace  is  a  satisfactory 
representation  of  the  field  trace,  or  if  the  parameters  should  be  adjusted  and  a  new 
model  trace  calculated.  This  process  is  repeated  until  a  maximum  number  of 
iterations  has  been  reached  or  the  model  fit  is  acceptable.  The  conventional 
methods  for  data  inversion  invert  all  parameters  simultaneously.  When  the  problem 
contains  a  large  number  of  parameters,  >diich  could  be  the  case  with  GPR  data,  the 
inv^ion  becomes  time  consuming  and  problems  with  convergence  arise.  A 
possible  inversion  technique  to  avoid  this  problem  is  to  use  a  layer  stripping 
method.  Layer  stripping  involves  inverting  for  parameters  one  layer  at  a  time  and 
then  removing  that  layer  and  proceeding  downward  through  the  subsurface.  The 
number  of  parameters  solved  for  at  a  given  time  are  reduced,  and  therefore  the  time 
required  to  perform  the  inversion  is  greatly  reduced.  The  layer  stripping  technique 
has  been  used  to  invert  dc  resistivity  data  (Pekeris  1940)  and  seismic  reflection  data 
(Yagle  and  Levy  1985). 


Manual  inversion  of  GPR  Data 


To  illustrate  the  procedure  involved  in  a  numerical  inversion,  the  manual 
inversion  of  GPR  field  data  is  presented  through  the  use  of  the  forward  modeling 
program  GPRMODv2  (Powers  and  Olhoeft  1995).  A  plot  of  the  field  data  is  given 
in  Figure  17.  The  hyperbolic  reflection  on  the  right  side  of  the  radar  section  is 
caused  by  a  4  feet  diameter  steel  pipe.  The  initial  parameter  values  chosen  to  invert 
the  data  are  listed  below.  Only  four  of  the  parameters  in  the  table  were  allowed  to 
vary;  all  others  were  assumed  constant.  The  thickness  of  the  material  overlaying  the 
pipe  was  calculated  based  on  an  assumed  velocity  of  0.  Im/s  and  a  one-way  travel 
time  of  18ns  (obtained  fi'om  plot  of  field  data). 


PARAMETER 

VARY 

VALUE 

Number  of  Layers 

no 

2 

Layer  Thickness 

yes 

h,=1.8m 

Conductivity 

no 

o,=10mS/m,  O2=10,000mS/m 

Eri 

yes 

Erii=8,  eri2=-l 

En. 

yes 

Ep.,=8, 

•Ce 

no 

t„=10ns,  t^=0ns 

no 

a„=0.8,  =1 

Pri 

no 

Pril“l>  PrI2~l 

Pt- 

no 

P^,  =  l,Pr_2=-l 

no 

x^,=0ns,  x^2=0ns 

“m 

no 

“h2=1 

Cr 

yes 

Offset 

no 

9.45  ns 
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Figure  18a  is  a  plot  showing  the  initial  model  fit  (solid  curve)  to  the  field  data 
(dotted  curve).  The  field  scan  is  the  GPR  trace  at  the  apex  of  the  pipe  reflection. 

The  GPR  field  data  (Figure  17)  show  that  the  onset  of  &e  pipe  reflection  occurs  at 
about  36ns  and  it  is  at  this  point  where  fitting  of  the  model  trace  will  begin.  The 
initial  model  requires  the  curve  to  be  expanded  vertically  and  shifted  to  the  left.  To 
accomplish  this,  the  dielectric  permittivity  is  increased  to  9.0,  giving  a  C  =1.5,  and 
the  l^er  thickness  decreased  to  1.34m  for  the  first  iteration  (Figure  18b).  The 
model  curve  is  now  in  approximately  the  correct  position  laterally,  however  the 
curve  still  requires  more  vertical  expansion.  For  the  second  iteration,  the 
permittivity  is  increased  to  10  with  C,=1.58  and  the  layer  thickness  must  be 
decreased  to  1.2m  (Figure  18c).  This  iteration  appears  to  have  smoothed  the 
variations  present  in  the  curve  between  55ns  and  70ns,  which  should  be  retained 
since  the  field  scan  shows  the  undulations  and  therefore  should  be  modeled.  In  the 
next  iteration,  the  permittivity  is  reduced  to  9.2,  C  =1.52,  and  the  layer  thickness 
increased  to  1.3m  (Figure  18d).  Now  the  model  curve  has  basically  the  same  shape 
as  the  field  scan,  so  the  dB  is  increased  to  7  to  vertically  enhance  the  model  response 
(Figure  18e).  Further  attempts  to  match  the  positive  peaks  in  the  model  curve  with 
Aose  in  the  field  scan  were  not  successful  so  the  model  fit  in  Figure  18e  is  accepted 
and  the  inversion  process  terminated.  The  parameters  used  in  each  iteration  step  are 
summarized  below. 


ITERATION 

-ErU 

JE^ 

-a. 

_h_ 

dB 

Initial 

8.0 

8.0 

1.41 

1.80 

0 

1 

9.0 

9.0 

1.50 

1.34 

0 

2 

10.0 

10.0 

1.52 

1.30 

0 

3 

9.2 

9.2 

1.52 

1.30 

0 

Final 

9.2 

9.2 

1.52 

1.30 

7 

Holding  many  of  the  parameters  constant  is  not  a  bad  practice  if  effects  caused  by 
dispersion  are  small  and  the  remaining  parameters  can  be  reliably  determined  or 
estimated.  Which  parameters  can  be  held  fixed  during  the  inversion  depends  on  the 
soil  properties  at  a  site.  The  iterative  procedure  is  considerably  simplified  when  the 
number  of  varying  parameters  involved  in  the  inversion  is  reduced. 


Chapter  5  GPR  Inverse  Modeling 


29 


6  Conclusions 


In  situations  where  ground  penetrating  radar  has  a  sufficient  depth  of 
investigation  capability  to  address  a  specific  problem,  it  has  the  greatest  horizontal 
and  vertical  resolution  of  any  geophysical  method  that  can  be  applied  to  the 
problem.  Existing  analysis  procedures  use  only  a  small  amount  of  the  information 
contained  in  the  received  GPR  waveform,  and  do  not  consider  electromagnetic  wave 
attenuation,  geometric  spreading,  dispersion,  scattering,  and  interface 
reflection/transmission  losses.  This  research  focused  on  the  assumptions  and 
considerations  necessary  for  developing  an  interpretation  procedure  that  utilizes  the 
full  GPR  waveform  for  determining  layer  thicknesses  and  the  EM  properties  of 
subsurface  materials. 

Forward  modeling  of  GPR  data  is  helpful  in  determining  if  it  is  feasible  to  use 
GPR  at  a  particular  site  given  certain  physical  properties  of  the  soil,  what  frequency 
antenna  is  required  to  detect  layers  of  varying  thickness,  and  what  type  of  response 
should  be  expected  for  various  combinations  of  conductivity,  dielectric  permittivity, 
and  magnetic  permeabihty.  The  forward  modeling  study  showed  that,  because  of 
the  superposition  of  primary  arrival  events  and  multiples,  arrival  times  other  than 
that  of  the  first  primary  reflection  cannot  always  be  determined  with  confidence  by 
visual  inspection  of  the  GPR  record.  The  application  of  data  processing  techniques 
to  detect  event  arrivals  that  are  too  closely  spaced  in  time  to  be  visually  resolved  is 
required. 

An  initial  investigation  of  the  inversion  of  GPR  data  was  undertaken.  There  are 
many  parameters  involved  in  the  inversion  of  frill  waveform  GPR  data.  The 
inversion  problem  would  be  unmanageable  without  making  certain  assumptions  and 
restraining  some  parameters  with  field  measurements.  Fortunately,  in  the  high 
frequency  range, lOOMHz,  dispersion  often  times  is  not  a  major  concern  so  the 
imaginary  component  of  the  dielectric  permittivity  and  magnetic  permeability  is 
zero.  If  the  soil  is  relatively  iron-free  so  the  magnetic  permeability  can  be  assigned 
a  fixed  value,  then  the  number  of  free  parameters  in  the  inversion  is  reduced  to  three 
per  layer,  the  real  dielectric  permittivity  in  GPRMODv2),  conductivity,  and 
layer  thickness,  plus  the  coupling  ratio  (C  =v'e/2). 

Both  crosscorrelation  and  deconvolution  were  shown  to  be  effective  at 
delineating  material  interfaces,  however  deconvolution  was  better  able  to  resolve 
interfaces  for  more  complex  transmitter  waveforms.  Identification  of  the  number  of 
interfaces  (hence  number  of  layers)  and  their  approximate  location  in  depth  and/or 
time  is  an  important  first  step  in  the  inversion  process. 
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RADAR  ANTENNAE 


Relative  Dielectric  Permittivity 


COLE-COLE  DISTRIBUTION 


Figure  2.  Cole-Cole  plots  of  the  real  and  imaginary  parts  of  the  dielectric 
permittivity  as  a  function  of  frequency 
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Figure  4.  Effect  of  angle  of  incidence  on  reflection  and  transmission  coefficients  . 
l^i  =  1^2  “  ^(b  propagation  constant  Yi  =  2,  Y2  =  4 


Figure  5.  Illustration  relating  antenna  separation  (s)  to  target  depth  (y).  Rule  of 
thumb  is  s=0.2y  (Annan  and  Cosway  1 992) 
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Figure  7.  Comparison  of 300,  500,  and  900  MHz  Ricker  waveforms  (left)  with 
calibrated  transmitter  (Tx)  waveforms  (right) 
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Figure  8.  Two-way  reflection  times  for  primary  and  selected  multiple  events  in  a 
concrete/water/sand  model  for  three  concrete  water^^ontents 
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Figure  9.  GPRM0Dv2  output  for  Case  172  (top)  and  Case  65  (bottom),  concrete 

w.c.  =  9.5%  and  900  MHZ,  wth  primary  reflections  (both  plots)  and  multiples 
(bottom  only)  identified 
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Figure  10.  GPRMODv2  output  for  Case  245  (top)  and  Case  137  (bottom);  concrete 

w.c,  =  0.0%  and  900  MHZ,  with  primary  reflections  (both  plots)  and  multiples 
(bottom  only)  identified 


Figure  11.  GPRM0Dv2  output  for  Case  149  (top)  and  Case  41  (bottom),  c'oncrete 

w.c.  =  9.5%  and  300  MHz,  with  primary  reflections  (both  plots)  and  multiples 
(bottom  only)  identified 
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Figure  12.  Comparison  of  concrete/air,  Cases  10-12,  with  concrete/water  reflections, 
Cases  28-30 
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Figure  14.  GPRMODV2  output  for  Case  252,  identifying  eight  Events:  two  primaty. 

reflections  and  six  multiples 
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Figure  1 8.  Illustration  of  manual  inversion  of  GPR  data,  a)  initial  model  fit  (solid 
curve)  to  field  data  (dotted  curve),  b)  model  fit  after  Iteration  1 , 
c)  model  fit  after  Iteration  2,  d)  model  fit  after  iteration  3,  e)  model  fit 
obtained  by  increasing  dB  level  ' 
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Figure  18.  (Continued) 
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Figure  18.  (Concluded) 
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Table  Al—Concrete/Air:  Cases  1-27 
Table  A2— ConcreteAVater;  Cases  28-36 
Table  A3-ConcreteAVater/Sand:  Cases  37-252 


Appendix  A  Tables  of  Cases  for  Forward  MexteTing  Probtem  Sets 


Table  A1 

Concrete/aIr:  Cases  1-27 


Table  A1  (Continued) 
Concrete/aIr:  Cases  1-27 


CM  Ico  1^  lio  |<o  ^  |« 


mm 


HiirnnH 

oirainio 


HflHin 


Table  A1  (Continued) 
Concrete/air;  Cases  1-27 
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Concrete/Water:  Cases  28-36 
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Table  A3 

ConcreteA/Vater/Sand  Variable  Water  Layer  Thickness  Test  Cases 
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Table  A3  (Continued) 

ConcreteAAfater/Sand:  Variable  Water  Layer  Thickness  Test  Cases  37-252 
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GPRM0Dv2  Forward  Modeling 
Results  for  Problem  Set  1 


Appendix  B  GPRMODv2  Forward  Modeling  Resuits  for  Problem  Set  1 
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Figure  B4.  GPR  model  calculation  results,  Cases  10-12 
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Figure  B9.  GPR  model  calculation  results,  Cases  25-27 
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Appendix  C  GPRMODv2  Fonward  Modeling  Results  for  ProWem  Set  2 
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Figure  Cl.  Problem  Set  2,  Cases  37-40,  multiples 


Figure  C3.  Problem  Set  2,  Cases  45-48,  multiples 
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Figure  C5.  Problem  Set  2,  Cases  53-56,  multiples 
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Figure  C9.  Problem  Set  2,  Cases  69-72,  multiples 
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Problem  Set  2,  Cases  73-76,  multiples 
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Figure  Cl  1.  Problem  Set  2,  Cases  77-80,  multiples 
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Figure  C13.  Problem  Set  2,  Cases  85*88,  multiples 
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Figure  C15.  Problem  Set  2,  Cases  93-96,  multiples 


Figure  C17.  Problem  Set  2,  Cases  101-104,  multiples 


Figure  Cl 9.  Problem  Set  2,  Cases  109-1 12,  multiples 


Figure  C21.  Problem  Set  2,  Cases  1 17-120,  multiples 


Problem  Set  2,  Cases  121-124,  multiples 


Figure  C23.  Problem  Set  2,  Cases  125-128,  multiples 


Figure  C26.  Problem  Set  2,  Cases  137-140,  multiples 


Figure  C27.  Problem  Set  2,  Cases  141-144,  multiples 
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Figure  C29.  Problem  Set  2,  Cases  149-152,  no  multiples 


Problem  Set  2,  Cases  153-156,  no  multiples 
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Problem  Set  2,  Cases  161-164,  multiples 
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Figure  C33.  Problem  Set  2,  Cases  165-168,  multiples 


Problem  Set  2,  Cases  169-172,  multiples 


Figure  C35.  Problem  Set  2,  Cases  173-176,  multiples 


Problem  Set  2,  Cases  177-180,  multiples 
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Figure  C37.  Problem  Set  2,  Cases  181-184,  multiples 


Problem  Set  2,  Cases  185-188,  multiples 
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Figure  C39.  Problem  Set  2,  Cases  189-192,  multiples 


Problem  Set  2,  Cases  193-196,  multiples 


Problem  Set  2,  Cases  201-204,  multiples 
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Figure  C44.  Problem  Set  2,  Cases  209-212,  multiples 


Mtef  Recompute  Shlft^Enter.Moln  M#nu 


Problem  Set  2,  Cases  217-220,  multiples 


Figure  C47.  Problem  Set  2,  Cases  221-224,  multiples 


igure  C48.  Problem  Set  2,  Cases  225-228,  multiples 


Problem  Set  2,  Cases  233-236,  multiples 


Problem  Set  2,  Cases  241-244,  multiples 


Figure  C53.  Problem  Set  2,  Cases  245-248,  multiples 


concr«te/Water/Sond:  w.c,-O.Ox 


Problem  Set  2,  Cases  249-252,  multiples 
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